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Abstract 

Dark energy models based on f(R) theory have been extensively 
studied in literature to realize the late time acceleration. In this pa- 
per, we have chosen a viable f(R) model and discussed its effects on 
the dynamical instability of expansionfree fluid evolution generating a 
central vacuum cavity. For this purpose, contracted Bianchi identities 
are obtained for both the usual matter as well as dark source. The 
term dark source is named to the higher order curvature corrections 
arising from f(R) gravity. The perturbation scheme is applied and 
different terms belonging to Newtonian and post Newtonian regimes 
are identified. It is found that instability range of expansionfree fluid 
on external boundary as well as on internal vacuum cavity is inde- 
pendent of adiabatic index T but depends upon the density profile, 
pressure anisotropy and f(R) model. 
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1 Introduction 



The "modified gravity" has become a standard terminology for the theories 
describing gravitational interactions which differ from the most conventional 
theory of general relativity (GR). In these modified theories, f(R) gravity is 
able to mimic the standard ACDM cosmological evolution and dark energy 
(DE) problem. Since the laws of gravity gets modified on large distances in 
f(R) models, this leaves several interesting observational signatures such as 
modification to the spectra of the galaxy clustering [lj, cosmic microwave 
background [2 J and weak lensing [3]. These models have some vacuum solu- 
tions with null scalar curvature that allow to recover certain GR solutions. 
In addition, it has many other applications such as inflation, local gravity 
constraints, cosmological perturbations and spherically symmetric solutions 
in weak and strong gravitational backgrounds. 

The most important feature of f(R) gravity is to provide the very natural 
gravitational alternative for DE without adding any matter component. Let 
us now show that how it can be related to the problem of DE by a straight- 
forward argument. When the Einstein-Hilbert (EH) gravitational action in 
GR, 

Seh = 7T I d A Xsf=gR, (1.1) 
2/c J 

is written in the modified form as follows 

S m = ^J d 4 x^f(R), (1.2) 

the addition of a non-linear function of the Ricci scalar demonstrates to cause 
acceleration for a wide variety of f(R) function, e.g., [I]- [13]. Variation of 
f(R) action with respect to the metric tensor leads to the following fourth 
order partial differential equations 

F{R)R aP -^f{R)g afi -V a VpF{R)+g a f } nF{R) = KT af s, {a, f3 = 0, 1, 2, 3), 

(1.3) 

where F(R) = df(R)/dR. Writing this equation in the form of Einstein 
tensor, it follows that [H] 

G a , = j(T^ + Tif), (1.4) 
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where 



rp(D) _ 1 



f(R) RF(R) ^ + VaVpF{R) _ gapUF{R) 



;i.s) 



Equation ( II .41) shows that effective stress-energy tensor plays the role 
of matter source in the field equations with purely geometrical origin. This 
approach may provide all the matter ingredients needed to tackle the dark 
side of the universe. Thus if we restrict "dark source" such that it does 
not satisfy the usual energy conditions then it can play the role of both 
dark matter and DE. Consequently, this theory may be used to explain the 
expansion of the universe and effects of DE in gravitational phenomena. This 
approach is sometimes convenient to use when we study the DE equation of 
state [15] as well as the equilibrium description of thermodynamics for the 
horizon entropy [16]. The trace of Eq. (ll.3j) determines the dynamics of the 
scalar field F(R) given by 

F(R)R~2f(R) + 3DF(R) = kT. (1.6) 

Notice that there exists a de Sitter point which corresponds to a vacuum 
solution (T = 0) at which the Ricci scalar is constant. Since QF(i2) = at 
this point, so we obtain from the above equation 

F(R)R-2f(R) = 0. (1.7) 

It is mentioned here that the model f(R) = SR 2 satisfies this condition and 
yields the exact de Sitter solution [T7]. In the model f(R) = R+SR 2 , because 
of the linear term in R, the inflationary expansion ends when 5R 2 becomes 
smaller than the linear term. This is followed by a reheating stage in which 
the oscillation of R leads to the gravitational particle production. It is also 
possible to use the de Sitter point for DE. 

The gravitational collapse is an important and long standing issue in GR. 
Recently, it has gained attention in modified theories as well. During gravi- 
tational collapse, self-gravitating objects may pass through phases of intense 
dynamical activities for which quasi-static approximation is not reliable. For 
instance, the collapse of very massive stars [18J, the quick collapse phase 
yielding neutron star formation [TJ] and the peculiar stars. The dynamical 
equations are used to observe the collapsing process while the vanishing ex- 
pansion scalar condition is developed in connection with the description of 
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voids. The vanishing of expansion scalar requires that the innermost shell of 
the fluid should be away from the center, initiating therefrom the formation 
of the cavity within the matter distribution [20j|2T]. This natural appearance 
of a vacuum cavity suggests that they might be relevant for the moddelling of 
cosmological voids. The Skripkin model is the first example satisfying expan- 
sionfree condition [22] • This model corresponds to evolution of spherically 
symmetric non-dissipating fluid distribution with constant energy density. It 
is also remarked that the expansionfree condition is a sufficient but not a 
necessary condition for the appearance of cavities. Cavities described under 
different kinematical condition are discussed by Herrera et al [23] . 

For the physically relevant models, the expansionfree evolution requires 
pressure anisotropy in the fluid distribution and inhomogeneity in the energy 
density. To assume isotropy of the pressure together with the expansion-free 
condition, we impose that the energy-density is independent on the time- 
like coordinate, which severely restricts the models [21]. A stellar model 
can exist only if it is stable against fluctuations. The problem of dynamical 
instability is closely related to the structure formation and evolution of self- 
gravitating objects. Chandrasekhar [22] was the first who worked in this 
direction. Afterwards, this issue has been investigated by many authors for 
adiabatic, non-adiabatic, anisotropic and shearing viscous fluids [26]-[30j. 
We have investigated the problem of DE and gravitational collapse in f(R) 
gravity [SUES]. 

In this paper, we are concerned with spherically symmetric stars having 
locally anisotropic fluid distribution inside and would investigate how f(R) 
terms affect the dynamical instability of expansionfree fluid evolution. The 
format of the paper is as follows. In section 2, basic equations are given. Sec- 
tion 3 is devoted to study the perturbation scheme and a well-known physical 
f(R) model. In section 4, Newtonian and post Newtonian approximations 
are taken into account and different terms belonging to these regimes are 
identified. Also, dynamical equations are investigated under the conditions 
of vanishing scalar and instability conditions of fluid evolution are discussed. 
The last section 5 provides the summary of the work. 

2 Field Equations and Dynamical Equations 

We consider a 3-D hypersurface E* 6 ), an external boundary of the collapsing 
spherically symmetric star, which divides a 4D spacetime into two regions 
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named as interior and exterior spacetimes. The interior spacetime to X* 6 ) is 
described by the most general spherically symmetric metric as follows 

dsl = A 2 (t, r)dt 2 - B 2 (t, r)dr 2 - C 2 (t, r)(d6 2 + sin 2 6d<p 2 ). (2.1) 

The exterior spacetime is described by the Schwarzschild metric given by 

ds 2 + = (l - dv 2 + Idrdv - r 2 (d9 2 + sin 2 9d(f) 2 ), (2.2) 

where M represents the total mass of the system inside the boundary surface 
and v is the retarded time. The fluid filling the spherically symmetric 
star is assumed to be locally anisotropic with inhomogeneous energy density. 
Thus the energy-momentum tensor for such a fluid is given by 

T al3 = (p + p±)u a up - p ± g a/3 + (p r - p±)xaXp, (2.3) 

where p is the energy density, p± the tangential pressure, p r the radial pres- 
sure, u a the four- velocity of the fluid and Xa is the unit four- vector along the 
radial direction. These quantities satisfy the relations 

u a u a = l, X a Xa = -l, X a u a = (2.4) 

and are obtained from the following definitions in co-moving coordinates 

u a = A" 1 ^, X a = B-'S?. (2.5) 

The expansion scalar, O, is given by 

e =«=i(! +2 §)< < 2e > 
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where dot and prime represent derivatives with respect to t and r respectively. 
For spherically symmetric interior metric, the field equations f ll.4j) become 



IB C\C 
~B + c)c 
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(2.7) 
(2.8) 



(2.9) 
A'B' 



AB 



F F" F 
A2 + B~2 + A* 



(2.10) 



The dynamical equations help to study the properties of collapsing pro- 
cess. To formulate these equations, the mass function introduced by Misner 
and Sharp is defined as follows [33] 



C 



C 



c 2 c 
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m(t, r) = -(1 + g^Cv) = _(! + _-_ 



(2.11) 



This equation provides the total energy inside a spherical body of radius "C" . 
From the continuity of the first and second differential forms, the matching of 
the adiabatic sphere to the Schwarzschild spacetime on the boundary surface, 
yields the following result j2Tj |3l] 



V(e) 

M = m(t,r) 



(2.12) 
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and 
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~2C 


( 2A 
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X 1 


-i) 


C 



(2.13) 



Using the field equations ( 12. 8 j) and (12. 9p in the above equation, we obtain 



Pr 



5 5 



+ 



01 



AB 



(2.14) 



As the expansion scalar describes the rate of change of small volumes 
of the fluid. It is interesting to note that expansionfree models define the 
two hypersurfaces, one separating the fluid distribution externally from the 
Schwarzschild vacuum solution while the other is boundary between the cen- 
tral Minkowskian cavity and the fluid. Taking £W [i stands for internal) 
to be the boundary surface of that vacuum cavity and matching it with 
Minkowski spacetime, we have 



m(t, 



0, -p, 



(D) 



11 



B 2 



+ 



T (D) 
J 01 

AB 



(2.15) 



The physical applications of such models lie at the core of astrophysical back- 
ground where the cavity within the fluid distribution is present. It may help 
to investigate the voids on cosmological scales [35]. Voids are the sponge- 
like structures and occupying 40%-50% volume of the universe. There exist 
different sizes of the voids, i.e., mini- voids [36] to super- voids [37]. Obser- 
vational data shows that voids are neither empty nor spherical. However, 
for the sake of investigations they are usually described as vacuum spherical 
cavities surrounding by the fluid. 

The proper time and radial derivatives are given by 



1__9 

Adi 



C 



l_d_ 



(2.16) 



where C is the areal radius of the spherical surface. The velocity of the 
collapsing fluid is defined by the proper time derivative of C, i.e., 



U = D T C = - 



(2.17) 
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which is always negative. Using this expression, Eq. (12.11j ) implies that 



E 



a_ 

~B 



u 2 + 



2m 
~C 



1/2 



(2.18) 



The rate of change of mass in Eq. (12.11]) with respect to proper time, with 
the use of Eqs. fl2Tj) -f l2~T0|) . is given by 



-K 



T 2F 



t (d) \ T {D) 



c 2 . 



(2.19) 



This represents variation of total energy inside a collapsing surface of ra- 
dius C. The presence of dark fluid components shows the contribution of 
DE having large negative pressure. These terms appear with positive sign 
representing negative effect, hence decrease the rate of change of mass with 
respect to time. Now, it depends upon the strength of DE terms that they 
may balance the positive effect of effective radial pressure or overcome on 
them. Likewise, we have 



Dcm 



K 

2F 



P + 



T (D) 

A 2 



UT, 



(D) 



01 



E AB 



C 2 



(2.20) 



This equation describes how energy density and curvature terms influence 
the mass between neighboring surfaces of radius C in the fluid distribution. 
Here the rate would decrease in the consecutive surfaces by the repulsive 
effect of DE. Integration of Eq.f l2.20p with respect to "C" leads to 



C ^2 

2F 



T, 



(D) 



P 



00 



UT, 



(D) 



01 



A 2 



E AB 



dC. 



(2.21) 



The dynamical equations can be obtained from the contracted Bianchi 
identities. Consider the following two equations 



(D) 



u, 



0. 



(D) 

rpaft _|_ rpa{3 



Xc 







(2.22) 



which yield respectively 



A + AB {P+P r / 
§ + + ^ 



+ ~(P + P±) + D 1 



C 

~BC 



= 0, 

^2 = 0, 



(2.23) 
(2.24) 
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where D\ and D 2 are components of dark source given in Appendix (15.11 



3 The f(R) Model and Perturbation Scheme 



Obtaining the correct dynamics of the background cosmological model is not 
sufficient for any theory to be viable. It is practically impossible to sepa- 
rate different f(R) models without using cosmological perturbations. In this 
section, we apply perturbation scheme on the field equation and dynamical 
equations in order to investigate the instability conditions of collapsing fluid 
evolution. We consider following f(R) model 



where 5 is any real number. When we take conformal transformation of f(R) 
action, it is found that the Schwarzschild solution is the only static spherically 
symmetric solution for the above model [38]. The stability criteria for this 
model is restricted to S > which corresponds to / (R) > 0. For 5 = 0, 
GR is recovered in which black holes are stable classically but not quantum 
mechanically due to Hawking radiations. Since such features also found in 
f{R) gravity, hence the classical stability condition for the Schwarzschild 
black hole can be expressed as / (R) > [39J. 

In our perturbation scheme, we assume that initially all the quantities 
have only radial dependence, i.e., fluid is in static equilibrium. After that 
all the quantities and the metric functions have time dependence in their 
perturbation. This is given by 



f(R) =R + 6R 2 , 




A(t,r 
B(t,r 
C(t,r 
p(t,r 
p r (t,r 
p±(t,r 
m(t, r 

e(t,r 









(3-2) 
(3-3) 
(3-4) 
(3-5) 
(3-6) 
(3-7) 
(3-8) 
(3-9) 
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Also, the Ricci scalar in f(R) model follow the same scheme as 



R(t,r) = R (r) + eT(t)e(r), 
f(R) = R (l + 25R ) + eT(t)e(r)(l + 25R ) 
F(R) = (l + 25R ) + 2e5T(t)e(r), 



(3.10) 
(3.11) 
(3.12) 



where < e < 1. By the freedom allowed in radial coordinates, we choose 
Co(r) = r as the Schwarzschild coordinate. 

The static configuration of the field equations fl2.7p - fl2.10p is obtained as 



1 + 25R 
1 

i 

1 + 2SR 
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i 
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A 

An 



5Rl 
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25R' 



Bl 
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o L^o 



A B r 



(A' 


B' \] 


Uo 


Bo)\ 



r A B 



(3.13) 



(3.14) 



(3.15) 



Applying the perturbed quantities in Eqs. fl3.2p -f l3.12p to the field equations 
f l2TP - fl2TIO|) along with Eqs.f l3TT3|) -f l3TT5|) . the resulting equations f l53|) -f l5T6|) 
are given in Appendix. In its static configuration, the dynamical equation 
(JOIl is identically satisfied while fl2T2^) yields 



A' 2 1 

PrO + (P0 + Pro)-T- + ~{PrO ~ P±o) + ~ 

Aq r k 



— 3-Rn 



35RoRn 



25R' B' 



An 



Bl Bo \A + r 



45R' ' 
TBI 



25R' Aq 



0. 



(3.16) 



The perturbed configurations of Eq. (12.231) leads to 



1 

Ao 



P + (Po + Pro)T^- + 2(p + Pl _ )f- + D 3 f 
Bo r 
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where D 3 is given in Appendix. Integrating this equation with respect to 
time, we get 



P 



b c 
(po + Pro)-^- + 2(p + Rlo)- + D 3 
B r 



T. 



(3.17) 



The perturbed part of Eq. (12.241) is obtained as 
bT 



Bl 



A' 2 

PrO + (PO + Pro)-^ + ~{PrO ~ P±o) 

A r 



(PO + Pro) 



a y T 
A~o) B~o 



+ 2(p r0 - p± Q ) (-) ^- + ^- + (p + Pr)^- + -(pr -P±) + D 4 = 
\rJ B B A r 

where a lengthy expression for D4 is written in Appendix. Perturbation on 
Eq.(J21ID yields 



m 



ni 



T 







Bl)' 




(<- 


b 


Bo 



(3.18) 



(3.19) 



From the matching condition Eq. fl2.14p with Eq.( l3.6[) . it follows that 



PrO 



Pr 



E w 5R 2 2SR' (A' q 2 



Bl \Ao r 



E w 25eT 25T 



Al A0B0 



-e 25R + 



+ 



3-Rq 



2(1 + 25.Ro) 

1 (K 2 



+ 



R'n fa' 2A' n 2c' 2c 46 



(3.20) 
Ro 



1 + 2SR 



+ 



+ 



+ 



Bl \A r 



\A A Bo r r 2 B r 
SeR'n 



1 + 2SR C 



(3.21) 



Substituting the above equations in Eq.f l5.5p . we obtain 



a(r)T(t) + P(r)T(t) + -y(r)T(t) = 



(3.22) 



where a, (3, and 7 are provided in Appendix. For the sake of instability re- 
gion, we assume that all the functions involved in the above equation are such 
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that a, (3 and 7 remain positive. The corresponding solution of Eq. ( I3.22|) is 
given by 



T(t) = — exp(a; S (e)£), where w E ( e ) = — - ~*~ - ^ ^ . (3.23) 

Here we impose that the system starts collapsing at t = —00 such that 
T(— 00) = 0, keeping the system in static position. Afterwards with the 
increase of t, it goes on collapsing. Using Eqs. fl3.9p and ( 12.61) . it follows that 

t / b 2c\ , 
6 = — — + — . (3.24) 



A \B r 

The expansionfree condition implies that 

be , 
— = -2-. 3.25 
B r 

It is interesting to mention that expansionfree fluid evolution causes a blowup 
of shear scalar which results the appearance of a naked singularity [21, 30J. 
The dynamical instability of collapsing fluids can be well discussed in term 
of adiabatic index T. We relate p and p r for the static spherically symmetric 
configuration by assuming an equation of state of Harrison- Wheeler type as 
follows pg SO] 

p r = r Pr0 p. (3.26) 

PO + PrO 

Here T measures the variation of pressure for a given variation of density. 
We take it constant throughout the fluid evolution. 



4 Newtonian and Post Newtonian Terms and 
Dynamical Instability 

This section provides help to identify the terms belonging to Newtonian (N), 
post Newtonian (pN) and post post Newtonian (ppN) regimes. This is done 
by converting relativistic units into c.g.s. units and expanding upto order 
c -4 in the dynamical equation. For the N approximation, we assume 

P0>Pr0, P0>P±0- (4.1) 
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For the metric coefficients expanded upto pN approximation, we take 



Gm 



~2 ' 



Bn 



1 + 



Gm 



-a ' 



(4.2) 



where G is the gravitational constant and c is the speed of light. From 
Eq.f l3.15p . we can write 
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+ rB 
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1 + 25R 
Also, from Eq. fl3.18p . we obtain 

K_ 

Bo t(t — 2m ) 
Using Eqs.f l3~T8|) and flBTTD . it follows that 

A' Q 2r 3 (K Pr0 -R - 361%) + A5r(Ro - 2R' r + iR' m ) + 4m 



(4.3) 



(4.4) 



An 



4r(r - 2m )(l + 25.Ro + SR' r) 



(4.5) 



Substituting Eqs. ( 14.2[) . ( 14.3p . (14 .5p in ( 13. 16}) and doing some algebra, the first 
dynamical equation in relativistic units is as follows 
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In view of dimensional analysis, this equation can be written in c.g.s. 
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units as follows 
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(4.7) 



Expanding upto c~ 4 order and rearranging lengthy calculations, the resulting 
expression given in Appendix. Here the order of c differentiates the terms 
in (I5.12p . (15.131) and (I5.14p respectively as follows 

terms of order c° correspond to N- approximation, (4-8) 
terms of order c~ 2 correspond to pN-approximation, (4-9) 
terms of order c -4 correspond to ppN-approximation. (4.10) 

This strategy may help to discuss some physical results in a certain regime 
by discarding terms belonging to other regimes. Applying the expansionfree 
condition (I3.25P on (I3.17p . we get 



P 



2(Pro-p±o)- + Da 
r 



T. 



(4.11) 



where D 5 is given in Appendix. This equation shows that the perturbed 
energy density depends on the static pressure anisotropy and higher order 
corrections. This fact supports the expansionfree condition that change in 
energy density depends exclusively on pressure anisotropy. Substituting this 



value of p in (I3.26p . we note that the expression p r and also p-^ in (14. 51) 
are of ppN order approximation. Thus we neglect them in the following 
calculations in order to discuss the instability conditions upto pN order. 

Applying expansionfree condition on the dynamical equation ( 15. 8ft along 
with value of T in (I3.23P and following the choice of radial functions 
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after tedious algebra, it follows that 
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Inserting Eqs. (14.21) . (14. 4D and (14.51) upto pN order (with c = G = 1) in the 
above equation, we obtain 



/ / \ 4co , w x ao 4<5eo 45eo^ 
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-(r - 4m )(a (r - ra ) + 8c r) 



H (r + 2m ) 



2r 



(r + 2mr 



x (2 + 3<H2 )] 1 + ^(r + 2m 



r 2 K 



25/?; 



-(r - 4m ) 



+ 



25R' 



(r — 2m r; 
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x [4c + a (r + 2m,)] + 2 -^^5R (r + 2m ) - 4 ^ 6 ° 

rot- 



K 



x (r - m )25R + 



% + 2m 2 ^r + 2m ' 



x [2e + 



m 2 ) + 12coi^] + 



32<5c 



r — 2m )[i?o(r — m )m 



45e c 



- 3m ) [r - 4c r 2 + m (r + 2m )] 
2r 3 (/tp r0 + 5i?o) + 46r(-2.R r + 4i? m ) + 4m " 



4r(r - 2m )(l + 25R + 5R' r) 



x 



ra (Po + Pro) H (?" + ™o) 



r 2^ 



V 



K 



(4c - a r) 



+ 



2<5e n 



ft 



2co 
rft 



x (r + 2m ) j^^(r-2m )l + ^-{r - 2m ){r + 2m ) + {a rR' 
I r J 1 /tr d 



+ e + 2<) + ^(r + 2m ) 2 {-10rc + - 2c (r + m ) 



r 5 ft 



a 



x (r + m )} - 2c Pro] ~ 



25u 2 



T ft 



(r + 2m )e 



4m^ 



+ (r + 2m ) 



x (e + 2c R' )] + %(r + m )e-<^>' + ^ [2c Qj R (r + 2m ) 



rft 



R' {r - 3m )(4c - a r) - e r + 



rft 



(r — m )(r + 2m ) 



x 



2r 3 (ftp r0 + £gg) + ASr(-2R' r + 4fl m ) + 4m 
4r (r - 2m ) (1 + 25.Ro + 5R' r) 



= 0. 



(4.14) 



In general, the instability range depends upon the index T as it measures 
the compressibility of the fluid. However, the above equation is independent 
of r which shows that instability region is totally depends upon the pres- 
sure anisotropy, energy density, chosen f(R) model and arbitrary constants. 
Notice that independence of T factor indicates that under expansionfree con- 
dition, fluid evolves without being compressed. In this way, the given f(R) 
model shows the consistency of the physical results with expansionfree con- 
dition. 

Notice that in the above equation, some terms appearing from ppN ap- 
proximation. In the following, we are interested in discussing the dynamic 
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instability of N-regime, so we ignore the terms belonging to pN and ppN 
approximations in Eq.f l4.14]) as 

i / i / \ 4co ao 4<5eo 
2co Pro + (Pro - P±o) I- (Po + Pro) (r + m ) 1 7 

+ ^^(r + 2m )-^(r + m )[(r + m )25R i 1 



x 



K,r° 
2cq 
rn 

25R' 



KT 



r 

4<5e 



r + 2m 
r — 4mo 



2r 



+ 2m )5R c 



.1 



+ ^-(r + 2m ) 



+ 



V - 2m )[4c + a (r + 2m )] 



r -m )(25i?o) + 



2p± 



r + m )e 



r — 3m )(l — 4cor) + 



i(U>: ' e) (r + 2m )(eo + 2c ^ 



7" ft 



2e R m 



r 3 K 



35R (r + 2m ). 



(4.15) 



In order to fulfill the instability of expansionfree fluids, we need to keep 
all the terms positive in Eq. (14.151) . Here, we assume that all the arbitrary 
constants and dynamical quantities are positive whereas p' r0 < showing 
that pressure decreases during collapsing process. In addition, we need to 
satisfy the following constraints 



Pro > Pud, — > c > 0, r > 4m . 
4r 



(4.16) 



Thus the system would be unstable in N-approximation as long as the above 
inequalities are satisfied because other constraints on m are followed by the 
last inequality in Eq. fl4.16p . 

For instance, if we assume that the scalar curvature is constant, i.e., 
Ro(r) = R c = constant and to = 0, then Eq. fl4.15p reduces to 

2c \p' r0 \ + (Pro - P_lo)— + ^^(r + 2m )(l + 2m' )5Ro 

r r z K 



+ (r + mo) 



(po+Pro)— + — (l+m' )8Ro+ 



(4.17) 
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Applying constant curvature condition on Eq. ( 12.21) ) and inserting in the 
above equation, we have 



o i ' i , ( \ Ac ° i RcC ° 

2c \p r0 \ + {PrO ~ P±0) H -^T~ 

r t z k 



R c + 



np r z 



1 + 25R C 
Por 2 dr 



(2 + 35R C ) + 



r + rR c + 
rR c 



[1 + 25 R c 

K 



E« 



Por 2 dr 



+ 



SCO 

2p ± 



, a 26R c c 

[PO+PrO) 1 

r nr 



2(1 + 26R C 

2 



R c + 



1 + 25.Ro 



-e s 



0. 



(4.18) 



Let us consider an energy density profile of the form po — Ar™, where A 
is a positive constant and — oo < n < oo. Substituting this value of po in 
Eq. jjUHD , it follows for n ^ -3 



O I'll/ N 4c , RcCQ 

2cq p r0 + (p r0 - p±o) h —5— 



/tAr 



n+2 



1 + 25R C 



(2 + 35R C ) + 



r + rR c + 
rR c 



3(1 + 25 R c ) 



r n+3 _ 



n+3 
EM 



X (r 
2p+ 



n+3 „"+3 



(Ar +p r0 ) 1 



0. 



2 ' 6(l + 2«? c ) 



1 + 2SR 



(4.19) 



Thus instability range depends upon the positivity of p r o — p+o and r 



n+3 



r n+3 

PrO > P±0 an d rn+3 > 

tonian regime, expansionfree collapse proceeds without compression as the 
adiabatic index does not involve in all the calculations. 



i.e., the system would hold the instability of expansionfree fluid for 

Finally, it is remarked that close to the New- 



„n+3 
E« 



5 Summary 

In this paper, we have studied the problem of gravitational collapse in f(R) 
theory which is strongly motivated by the observational data collected from 
supernova. Here the higher order curvature terms are thought to be the origin 
of DE causing acceleration which are treated as the matter part of the field 
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equations. We have considered the spherically symmetric collapsing stars 
made up of locally anisotropic fluid and evolving under the expansionfree 
condition. The curvature terms appear to affect the passive gravitational 
mass and rate of collapse. Also, f(R) DE slows down the rate of collapse 
due to its repulsive effect. 

The dynamical equations help to investigate the evolution of gravitational 
collapse with time and yield the variation of total energy inside a collapsing 
body with respect to time and adjacent surfaces. We have formulated these 
equations by using contracted Bianchi identities both for the usual mat- 
ter and effective energy-momentum tensor independently for a well-known 
f(R) = R + 5R 2 model. The first dynamical equation is used to identify the 
terms belonging to Newtonian, post Newtonian and post post Newtonian 
regimes. We have used the concept of relativistic and c.g.s units. The second 
dynamical equation is used to discuss the instability range of expansionfree 
fluid evolution upto pN order. 

Perturbation scheme is applied on the field equations and dynamical equa- 
tions. The study of resulting equations shows that the range of instability 
is independent of adiabatic index T which generally plays central role in the 
definition of instability range. For example, for a Newtonian perfect fluid, 
the system is unstable for T < 4/3. In our results, the independence of T 
shows the consistency of expansionfree condition with f(R) gravity because 
this condition requires that fluid would evolve without compressibility. More- 
over, the instability range depends upon the anisotropy of radial pressure, 
energy density and some constraints which arise for keeping the positivity of 
the dynamical equation in Newtonian approximation. Assumption of con- 
stant scalar curvature implies that the above dependence of instability range 
also describes the instability of cavity itself with the additional information 
that r n+3 > r^Jlf should be satisfied. It is mentioned here that at pN regime 
only relativistic effects are taken into account, however, physical behavior of 
the dynamical equation would be the same. 
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Perturbed field equations: 
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